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General relativistic simulations of black hole-neutron star mergers have currently been limited to low-mass 
black holes (AIbh < 7Mq), even though population synthesis models indicate that a majority of mergers might 
involve more massive black holes (AIbh > IOMq). We present the first general relativistic simulations of black 
hole-neutron star mergers with A/bh ~ 10A4q. For massive black holes, the tidal forces acting on the neutron 
star are usually too weak to disrupt the star before it reaches the innermost stable circular orbit of the black 
hole. Varying the spin of the black hole in the range aBu/MBu = 0.5 - 0.9, we find that mergers result in the 
disruption of the star and the formation of a massive accretion disk only for large spins aBH/AfsH > 0.7 - 0.9. 
From these results, we obtain updated constraints on the ability of BHNS mergers to be the progenitors of short 
gamma-ray bursts as a function of the mass and spin of the black hole. We also discuss the dependence of the 
gravitational wave signal on the black hole parameters, and provide waveforms and spectra from simulations 
beginning 7-8 orbits before merger. 

PACS numbers: 04.25.dg, 04.40.Dg, 04.30.-w, 47.75.+f, 95.30.Sf 



I. INTRODUCTION 

The study of black hole-neutron star (BHNS) mergers of- 
fers an opportunity to observe a general relativistic system in 
the strong field regime, in the presence of material at supernu- 
clear density, magnetic fields, shocks, and intense neutrino ra- 
diation. As for other compact binaries, recent interest in black 
hole-neutron star systems is due in large part to their potential 
as sources of gravitational waves detectable by the next gen- 
eration of ground-based detectors. Advanced LIGO IH] and 
VIRGO Indeed, whether for detection or parameter esti- 
mates, these detectors require theoretical templates in order to 
extract the gravitational wave signal from the detector noise. 
And in the orbits immediately preceding the merger, as during 
the merger itself, these binaries can only be properly modeled 
by numerical simulations in general relativity. 

The emission of gravitational waves is not the only interest- 
ing feature of BHNS mergers, however For compact binaries 
containing at least one neutron star the outcome of the merger 
is by itself an important question. The formation around the 
remnant black hole of a massive, thick, and hot accretion 
disk from the remains of a tidally disrupted neutron star of- 
fers ideal conditions to power short-hard gamma ray bursts 
(SGRB). This possibility, combined with the association of 
SGRBs with relatively old stellar populations, makes BHNS 
and NS-NS mergers some of the most likely progenitors of 
SGRBs (see e.g. yj] for a review). 

Black hole-neutron star binaries have not yet been as widely 
studied as NS-NS or BH-BH systems. The first simulation 
of a BHNS merger in general relativity was done in 2006 by 
Shibata & Uryu Since then, various groups have studied 
the influence on the dynamics of the merger and the gravi- 
tational wave signal of the black hole spin iH-Ql, the mass 
ratio ins, 01 J the neutron star equation of state ifTl- floll . the 
eccentricity of the orbit iflTIl , and the presence of a magnetic 
field llT2ll . A review of those results can be found in ifisll . 
These simulations indicate that mergers of low-eccentricity 



BHNS binaries show two potential qualitative behaviors. In 
the first, the neutron star is tidally disrupted before it reaches 
the innermost circular orbit (ISCO) of the black hole. Ma- 
terial from the star is then either quickly accreted onto the 
black hole, or ejected in a long tidal tail. As material in the 
tail falls back onto the black hole, it then forms an accretion 
disk. These disks are usually fairly thick, have temperatures 
of a few MeV, can contain a significant fraction of the initial 
mass of the star (a few tenths of a solar mass), and maintain a 
baryon-free region along the rotation axis of the black hole in 
which relativistic jets could be launched. They therefore offer 
promising conditions for the generation of SGRBs. The sec- 
ond scenario lets the star reach the ISCO before it overflows 
its Roche lobe. No tidal tail is formed, and the black hole and 
neutron star merge directly. Numerical simulations have con- 
firmed that disk formation is favored by a high black hole spin 
(as the ISCO is closer to the black hole), a large neutron star 
(the star is more easily disrupted), and a small black hole mass 
(the size of the Roche lobe at the ISCO relative to the radius 
of the neutron star is smaller for lower mass black holes, thus 
making Roche-lobe overflow easier). 

All existing simulations of BHNS mergers in general rela- 
tivity have studied black holes of mass between 2 and 5 times 
the mass of the neutron star For such configurations, disk for- 
mation is the most likely scenario. It can only be prevented for 
anti-aligned spins Isl Ml, and compact stars or massive non- 
spinning black holes jsi [Toll . However, the choice of such 
low mass ratios is not very well motivated. Population syn- 
thesis models tend in fact to favor distiibutions in which a 
large fraction of BHNS binaries have a more massive black 
hole (AfsH > IOMq), especially for binaries formed in low- 
metallicity environments ifTil [Tsll . These predictions come 
with large uncertainties, but it still appears important to study 
the behavior of BHNS systems for higher mass black holes. 

Simulations of binary systems with high mass ratios are 
typically more costly than their equal-mass counterparts. The 
reason is that before disruption the time step of the evolu- 
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tion is limited by the minimum spacing of the numerical grid 
(Courant-Friedrichs-Lewy condition), which scales as the size 
of the smaller object. High mass ratio simulations thus require 
significantly more time steps per orbit if we want to maintain 
the same accuracy. Mergers are also more difficult in this case: 
the rapid accretion of matter onto the black hole from a rela- 
tively small object can only be resolved if the numerical grid 
around the black hole has an accordingly small spacing, and 
a small time step is thus required as well. These increasing 
costs, plus the fact that lower mass ratios offer richer physical 
effects, have led to a focus on low black hole masses in general 
relativistic simulations of BHNS mergers up till now. High- 
mass ratio BHNS binaries have, however, been evolved using 
approximate treatments of gravity (e.g. in lfl6l - [l9ll V These 
simulations tend to indicate that high spins are required for 
disk formation to be possible. However, approximate simula- 
tions have been known to disagree with results in full general 
relativity in the past, especially on predictions regarding the 
mass of the accretion disk or the amount of unbound material 
ejected in the tidal tail. 

In this paper, we focus on the features of BHNS mergers 
for black holes of mass 7 — IOA/0. In Section we re- 
view our numerical code and the specific challenges related 
to the evolution of high mass ratio systems, and describe re- 
cent improvements to our code. We also present updated di- 
agnostics regarding the accuracy of our gravitational wave- 
forms, and the influence of gauge choices on our results. Sec- 
tion HU] describes the initial conditions for our simulations, 
while Section HV] presents our numerical results. In particu- 
lar, we show that for massive black holes (i\/BH > IOMq), 
disks cannot be formed unless the spin of the black hole is 
high {aBH/MBa > 0.7). We also discuss the imprint of 
the binary parameters on the gravitational wave signal. Fi- 
nally, Section IV] focuses on consequences for the observation 
of SGRBs. 



II. METHODS AND DIAGNOSTICS 

The black hole-neutron star binaries presented in this pa- 
per are evolved using the numerical code developed by the 
SpEC collaboration [20]. The coupled system formed by Ein- 
stein's equations and the general relativistic hydrodynamics 
equations is solved using the two-grid method described in 
Duez et al. ll2lll . with the improvements discussed in Foucart 
et al. We refer the reader to ll^l for a detailed descrip- 
tion of the numerical methods used in SpEC, and only discuss 
here some recent modifications improving the general perfor- 
mance of the code, and facilitating the merger of high mass 
ratio BHNS binaries. 

In the two-grid method, Einstein's equations are solved 
within the generalized harmonic formalism using pseu- 
dospectral methods, while the fluid equations are evolved on 
a separate finite difference grid. This allows us to take advan- 
tage of the efficiency of pseudospectral methods for the evo- 
lution of Einstein's equations in regions in which the solution 
is smooth (i.e., away from the neutron star), while limiting the 
extent of the finite difference grid to regions in which matter 



is present. The price to pay is that, as the two sets of equations 
are coupled, source terms have to be interpolated between the 
two grids at each time step. 

The relativistic hydrodynamics equations are solved in con- 
servative form: the evolved conservative variables U satisfy 
equations of the form 

dtU + W-FiU) = S{U), (1) 

where the fluxes F and source terms S are functions of the 
variables U but not of their derivatives. A conservative shock- 
capturing scheme is characterized by a reconstruction method, 
whereby fluid quantities on either side of cell faces are con- 
structed, and a Riemann solver that supplies the resulting 
fluxes across each cell face. For these simulations, we use 
5th-order WENO reconstruction ll23H24ll with the smoothness 
indicator proposed by Borges et al ESll and HLLE fluxes ll26ll . 

Since the publication of Foucart et al. |@], other modifica- 
tions of the SpEC code have significantly improved the effi- 
ciency of our simulations. In the following sections, we dis- 
cuss the most important of these recent changes. 

A. Adaptive time stepper 

In previously published simulations, black hole-neutron 
star binaries were evolved using a fixed time step At = 
CAxinin, where Axmin is the smallest grid spacing and 
C some constant chosen so that the evolution satisfies 
the Courant-Friedrichs-Lewy (CFL) stability condition at all 
times. As the truncation error in these simulations is usually 
dominated by the effects of the spatial discretization, this is 
sufficient for the effects of the discretization in time to be neg- 
ligible. However, the resulting time step can be significantly 
smaller than necessary for a large fraction of the evolution. In- 
deed, the maximum value of C that gives stable evolutions can 
change over time' and is a priori unknown for the evolution 
of Einstein's equations on a pseudospectral grid. Additionally, 
the only way to confirm that the error in a given simulation is 
indeed dominated by spatial truncation error is to rerun it with 
a different choice for C. Using an adaptive time stepper, we 
instead choose at any time the largest time step for which the 
time stepping error is smaller than some tolerance. If the CFL 
instability begins to grow, this shows up as an increasing time 
stepping error and the time stepper automatically reduces the 
time step. In practice, time step control is done by evaluating 
the time stepping error through comparisons of the result of 
the time evolution with what would be obtained if a lower or- 
der scheme were used. We evolve using the 3rd-order Runge- 
Kutta algorithm, and compare with the results of a second- 
order method. The time step hn+i is then determined using a 
proportional-integral (PI) stepsize control ll27l - [30ll : 

(2) 



The variation in time of the CFL stability condition occurs mainly because 
the pseudospectral grid follows the evolution of the binary and, in particu- 
lar, contracts as the binary spirals in. 
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where e„ is the error measured at time step n, S = 0.9 is a 
safety factor, P « 0.13 and a « 0.23. 

For the evolution of a BHNS binary over 7-8 orbits, choos- 
ing the optimal At at all times can reduce the simulation time 
by a factor of ^ 1.5 — 2. Additionally, as the minimum value 
of C giving stable evolutions is unknown, our choice when 
using fixed time steps was usually overly cautious. The actual 
gain in the time required to evolve BHNS systems with SpEC 
at a given resolution is thus closer to a factor of 2 - 3. 

B. Control system 

The evolution of a black hole on a pseudospectral grid re- 
quires the excision of the region surrounding the singularity. 
As a consequence, the numerical domain has an inner bound- 
ary (Sin, typically a sphere in the coordinates of the numerical 
grid. If any information can enter the domain through Sin, 
an unknown boundary condition has to be imposed there. In- 
side the apparent horizon of the black hole, this can however 
be avoided for 'good' choices of iSin. More specifically, we 
compute the local characteristic speeds of the generalized har- 
monic equations (which are hyperbolic), and require that on 
Sin all characteristic speeds point out of the numerical grid. 
We also require that iSi„ remain inside the apparent horizon of 
the black hole. In practice, this can be difficult when the black 
hole is rapidly evolving through tidal distortion or matter ac- 
cretion. During merger, we try to keep the apparent horizon 
nearly spherical in the coordinates of the numerical grid, and 
slightly outside Sin, through the use of a coordinate map be- 
tween the grid coordinates and the inertial frame: 

r = ril + J2Yi„,ie,cl))cim), (3) 

Ini 

where r is the distance from the black hole center in the iner- 
tial frame, r the same distance in the grid coordinates, {0,(j>) 
the usual angles in spherical coordinates (in both frames), and 
cim arbitrary coefficients used to control the location of Sm 
in the grid frame. The quantity cqq thus controls the size of 
Sin, while the other coefficients control its shape. In previous 
simulations, the control system attempted to choose the cim 
so that the apparent horizon remained a sphere of fixed radius 
in the grid frame. This proved inadequate for high mass ratio 
simulations: there is a priori no way to know at what distance 
the apparent horizon should be kept from Sin, and making the 
wrong choice can easily lead to characteristic speeds entering 
the domain on Sin- Instead, we now choose cqo by using the 
characteristic speeds themselves as input for the control sys- 
tem. Indeed, dtCQo affects the velocity of the surface Sin in 
inertial coordinates, and can thus be used to increase or de- 
crease the characteristic speeds at will. The other coefficients 
are still chosen so that the apparent horizon remains spherical. 

As the time scale over which Sin needs to change varies by 
orders of magnitude between the early orbits and the merger, 
the control system used to fix the c;„i has also been modified 
in order to choose adaptively the time scale over which it at- 
tempts to damp any deviation of the motion of Sin from its 
desired behavior 



C. Adaptive Mesh Refinement on the spectral grid 

Another complication arising for high mass ratio mergers is 
the rapid motion of the high-density neutron star core across 
the numerical grid as it falls into the black hole. This high- 
density material generates large source terms in Einstein's 
equations, which vary quickly and can have very steep gra- 
dients. Maintaining a grid resolution high enough to resolve 
these features in all parts of the numerical domain where they 
might occur, and during the entire simulation, is far too costly 
to be realistic. In practice, we need to adapt the numerical 
domain at regular time intervals in order to increase resolu- 
tion in regions where the neutron star core is located. This 
was already done on the finite difference grid, by limiting the 
grid to regions in which matter is present i^. The possibil- 
ity of modifying the pseudospectral grid during a simulation, 
on the other hand, is a new feature of our code which was 
initially developed to improve the accuracy of black hole bi- 
nary mergers. The choice of how to modify the resolution 
in any given subdomain of the numerical grid relies on esti- 
mates of the truncation error based on the fall-off rate of the 
coefficients of the spectral expansion. For a smooth solution, 
these coefficients should fall off exponentially. In BHNS sim- 
ulations, the solution is not smooth, but the metric quantities 
are still continuous, and their spectral expansion converges as 
a power law. The determination of the required resolution is 
less reliable than in vacuum — and better methods to control 
that choice should be found if we want to take full advantage 
of the method — but in its current form this spectral adaptive 
mesh refinement (AMR) technique has already allowed us to 
evolve configurations that were otherwise too costly to sim- 
ulate accurately. In this paper, spectral AMR is used for the 
merger of case Q7S5, the only configuration in which the neu- 
tron star merges without any significant tidal disruption. 



D. Gravitational Wave Accuracy 

Gravitational wave templates for ground based detectors 
such as Advanced LIGO require accurate waveforms. Even 
the longest and most accurate simulations of binary black 
holes combined with Post-Newtonian results to generate hy- 
brid waveforms might currently be insufficient for precise pa- 
rameter estimates (i.e., the accuracy of the parameters might 
be limited by the qualit y of the templates instead of the noise 
in the LIGO data) ll3ll l32ll . Templates could well require a 
phase accuracy in the numerical waveform of 0. 1 rad over sig- 
nificantly more than the 30 cycles of the longest available nu- 
merical waveform ll32ll . Accuracy requirements are lower for 
detection purposes, and some parameters could be usefully 
constrained with waveforms at the current level of accuracy 
(e.g., the neutron star radius from NS-NS mergers ^^). But 
it is clear that template accuracy is a significant concern for 
gravitational wave astronomy. 

Binary neutron stars and BHNS binaries are particularly 
challenging, as the need to evolve the neutron star fluid adds 
new sources of errors to the simulations. In previous pa- 
pers Sl^, we evolved the system for only 2-3 orbits before 
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FIG. 1 : Phase error Sff) for the (2,2) mode of the gravitational wave 
signal as a function of the number of orbits, for simulation Q7S7. 
The wave is extracted at finite radius R = lOOMtotai, and we do not 
attempt to shift it in phase or time. 
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FIG. 2: (2,2) mode of the strain h for simulation Q7S7 at our two 
highest resolutions, extracted at finite radius R = lOOAftotai. We 
shift the signal in time and phase so that the phase at the peak of the 
signal is = 0, and the time t = 0. 
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merger, and the phase error in the gravitational wave signal 
was of the order of a few tenths of a radian even before merger. 
The simulations presented here are significantly longer (7-8 
orbits), and more accurate (the 'low' resolution used here is 
roughly equivalent to our standard resolution in J^). We thus 
need to reassess the accuracy of our waveforms. 

The error in the gravitational wave signal is evaluated using 
the Q7S7 simulation {q = 7, aBH/AfsH = 0.7), which we 
run at 3 different resolutions. The gravitational wave signal 
is extracted at finite radius R = 1007\/tot, where Mtot = 
Mbh + A/ns- We directly compute the strain h expended in 
spherical harmonics {h — y^, himYim), using the Regge- 
Wheeler-Zerilli method ll34l - i37ll . In Fig.[T] we show the raw 
phase difference between our high resolution simulation and 
the lower resolution simulations for the dominant /i22 mode. 
The phase is directly extracted from the real and imaginary 
parts of /i22, without any time or phase shift. The phase error 
in the low resolution case reaches 0.2 rad after only 2.5 orbits, 
as expected from the fact that it is equivalent to the resolution 
used in i^. At higher resolution, however, the same phase 
difference is only reached after approximately 4.5 orbits, and 
the overall accuracy is significantly better up to merger even 
though by then the phase shift is of the order of a radian. 

Such errors are typical of current general relativistic simu- 
lations. We can, for example, compare our results with those 
of the most recent BHNS mergers of Kyutoku et al. [7] . There, 
the error in the waveform is measured after applying a shift in 
time and phase. This is justified from an observational point 
of view by the fact that a detector would not be sensitive to 
such shifts, even though it might also give an optimistic view 
of the uncertainties at the time of merger Fig. |2] shows the 
wave from the two highest resolution of simulation Q7S7, but 
shifted so that the time at which the signal peaks is < = and 
the phase at that time is = 0. The errors appear comparable 
to Fig. 25 of 17], and significantly smaller than the differences 
due to changes in the parameters of the binary, discussed in 
Sec. UVBl 



E. Gauge dependence 

In the generalized harmonic formulation, the gauge is spec- 
ified through a choice of the variable Ha ~ gab^ c^'^x^ ■ For 
most runs reported here, we evolve this function as in our pre- 
vious work |6]. During the inspiral, we fix Ha in the comov- 
ing frame at its initial value. During the merger, we continue 
to fix Ha near the black hole; away from the hole, we damp it 
to zero. Because this gauge is presumably not optimal for the 
merger evolution near the black hole, significant grid stretch- 
ing is expected and found during this part of the evolution. To 
test whether the gauge behavior significantly affects the accu- 
racy of the simulations, we reran the merger phase of the Q5S5 
runs using the damped-wave gauge introduced by Szilagyi, 
Lindblom, and Scheel Issll . The resulting metric evolution 
is, of course, different (and noticeably more moderate), but 
the gravitational wave signal should be the same. To within 
numerical errors, this is indeed what we find. Also, the differ- 
ence in waveform between resolutions is comparable for each 
gauge, indicating that grid distortion is not the major source 
of waveform error 



III. INITIAL CONFIGURATIONS 

We generate initial data for our simulations with the elliptic 
solver SPELLS |f39l]. To obtain realistic initial configurations 
for the evolution of black hole-neutron star binaries, we have 
to satisfy the constraints in Einstein's equations, as well as 
choose a realistic initial state for the fluid forming the neu- 
tron star For the fluid we require hydrostatic equilibrium as 
well as an irrotational velocity profile. The procedure for solv- 
ing for these conditions and obtaining compact objects with 
the desired masses, spins and orbital parameters is detailed in 
Foucart et al. ll40ll . The constraint equations are solved in the 
extended conformal thin sandwich formalism ll4lll42ll . In that 
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formalism, the metric is decomposed as 

ds^ = g^^dx>'dx'' ~a'^dt'^+(j>%j{l3'dt+dx'){l3^dt+dx^) 

(4) 

where a is the lapse function, /3* the shift vector, 7^ the con- 
formal 3-metric and (j> the conformal factor The constraints 
can then be expressed as elliptic equations for (acj), f3'^ , (/)). 
The choice of the conformal metric, the trace of the extrin- 
sic curvature K^^^ — ~^Cg^y, and their time derivatives is 
then arbitrary. Together with the boundary conditions and the 
stress-energy tensor, these arbitrary choices will determine the 
characteristics of the spacetime under consideration. A stan- 
dard choice in the literature is gij = Sij, K = 0. However, 
for high black hole spins this leads to configurations far from 
equilibrium, and potentially large modifications of the black 
hole spin and mass at the beginning of the simulation. We pre- 
fer to choose the free data (gij , K) according to the prescrip- 
tion of Lovelace et al. ll43ll : close to the horizon of the BH, the 
metric matches the Kerr-Schild solution. (For the adaptation 
of that method to BHNS systems, see llioll .) 

In this paper, we are interested in the behavior of BHNS bi- 
naries for high mass ratio q = A/bh/^^ns = 7. We consider 
4 different configurations, summarized in Table U The first is 
at a slightly lower mass ratio g = 5, a range of masses already 
explored by Etienne et al. lUt], Kyutoku et al. ^ and Chawla 
et al. 1121 . The spin of the black hole is aBH/AfsH — 0.5, 
aligned with the orbital angular momentum of the binary. The 
other three simulations use g = 7. We show in Sec. lIVAl that 
these cases span the range of BH spins over which the quali- 
tative features of the merger vary: oehZ-^bh = 0.5, 0.7, 0.9. 
All configurations are started approximately 8 orbits before 
disruption, and have low eccentricity e < 0.005 (obtained us- 
ing the iterative method developed by Pfeiffer et al. lliill ). For 
the equation of state of nuclear matter, we choose a polytrope: 
the pressure P and internal energy e are expressed as functions 
of the baryon density po and temperature T\ 

P = npl+p^T (5) 



with r = 2 and k chosen so that the compaction of the star 
is C = M/R = 0.144. This corresponds to a stellar radius 
R ^ 14.4 km for a 1.4Mq neutron star Recent predictions 
regarding the radius of neutron stars indicate that this is prob- 
ably close to an upper bound on the radius of real stars (see, 
e.g., H). 

We use three different resolutions: 80^, 100^ and 120^ 
points on the finite difference grid (which only covers the 
neutron star), and 56.9"^, 64.5'^ and 72.0"^ points on the pseu- 
dospectral grid. At the highest resolution on a moderate 
size cluster, 8 orbits require approximately 40, 000 CPU-hrs 
(about a month on 48 processors). During mergers, the reso- 
lution of both the finite difference and pseudospectral grids is 
higher, and depends on the configuration studied. 



TABLE I: Initial configurations studied. aBn/AfBH is the dimen- 
sionless BH spin, A/bh.ns the ADM masses of equivalent isolated 
BH and NS, -Rns the radius of the star, O^^at the initial angular fre- 
quency of the orbit, Mtot the ADM mass of the system at infinite sep- 
aration, and e the initial eccentricity. The number of orbits A^orbits 
is measured at the point at which 10% of the mass has been accreted 
by the black hole. 



Name 




^BH 


Q _ -'^''nS 

-Rns 




Aorbits 


e 


Q7S5 


7 


0.5 


0.144 


3.84e-2 


1.1 


0.001 


Q7S7 


7 


0.7 


0.144 


4.14e-2 


8.3 


0.004 


Q7S9 


7 


0.9 


0.144 


4.47e-2 


8.8 


0.005 


Q5S5 


5 


0.5 


0.144 


3.55e-2 


7.2 


0.003 



IV. RESULTS 
A. Mergers 

Most existing simulations of BHNS mergers with lower 
mass black holes (A/bh 3 - 7Mq) show very similar qual- 
itative behavior As the neutron star gets close to the black 
hole, tidal forces cause the star to overflow its Roche lobe. 
Most of the matter is either rapidly accreted onto the black 
hole, or ejected into an extended tidal tail. Material in the tidal 
tail then rapidly forms a thick and hot accretion disk around 
the black hole (T ^ 1 - 5 McV). To obtain different behav- 
iors, one must consider either highly eccentric orbits iflTll — in 
which case partial disruption of the neutron star is possible — 
or anti-aligned black hole spins lUt] , where it becomes possible 
to directly accrete the entire neutron star since it reaches the 
innermost stable circular orbit before overflowing its Roche 
lobe. Negligible disk masses are also observed for nonspin- 
ning black holes and fairly compact stars at the high end of 
the range of studied BH mass (g = 4 - 5) isllTtl. 

The situation changes significantly when considering more 
massive black holes, g ~ 7. For such binaries, the standard 
behavior appears to be the direct merger of the two compact 
objects, without significant tidal disruption of the stars. The 
simulations presented in this paper use neutron stars with radii 
on the high end of the theoretically acceptable values. This is 
known to favor tidal disruption ^}^, as material on the sur- 
face is less tightly bound to the neutron star Even so, in our 
lower spin simulation {q = 7, aBH/AfsH = 0.5), 99% of 
the material falls into the black hole in less than 2 ms. Snap- 
shots of the matter configuration for that simulation are given 
in Fig. [3] Once accretion begins, most of the neutron star ma- 
terial remains within a coordinate distance of less than about 
40 km from the center of the black hole, and only a negligible 
amount of matter forms a tidal tail (approximately 0.5% of 
the NS mass). As spin misalignment, a lower spin and a more 
compact neutron star all work against tidal disruption, this 
gives us a strict lower bound on the black hole spin required 
to disrupt the star for black holes with mass A/bh > 7il/Ns- 

When the spin of the black hole is increased to 
obh/A/bh = 0.7, the disruption of the neutron star becomes 
clearly visible (Fig.|4] left). A long tidal tail extending more 
than 100 km away from the black hole and containing about 
8% of the initial mass of the star is created from the ejected 
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FIG. 3: Matter distribution for simulation Q7S5. Left: Most of the neutron star material plunges directly into the black hole at the time of 
merger. Density contours: po ~ (10^^, 10^'', 10^^, 10^^, 10^*^) g/cm? . Right: 1ms later, only 2% of the neutron star material remains outside 
the black hole, and most of it will rapidly accrete onto the hole. Density contours; po = (10^^, 10^^, 10^") g/cra? . 





material. About 5 ms after disruption, a proto-accretion disk 
forms (Fig. |4] center). However, the disk contains only a 
fairly small portion of the ejected material, and is still thin- 
ner {H/R ~ 0.1, where H is the half-thickness of the disk) 
and colder (T < 1 McV) than what is usually observed in 
mergers with less massive black holes. 

Rapid accretion from the tidal tail prevents the disk from 
reaching an equilibrium profile for about 20 ms. By that point, 
only 2.5% of the baryonic mass remains outside the black hole 
(Fig. m right). As the radius of the disk is about 100 km, its 
density is much lower than in the case of higher black hole 
spins (see below) or lower black hole masses: the maximum 
density is only about 2 x 10^° g/cm'^. The thickness {H/R ^ 
0.3) and temperature (T ^ 2 MeV) of the disk are closer to the 
equilibrium values found for g = 3 |@], but this still appears 
to be a less desirable configuration for the potential generation 
of a short gamma-ray burst. 

This is no longer the case for the high spin configuration, 
obh/A^bh = 0.9. The disruption (Fig.|5] left) is very similar 
to the previous case, with a long tidal tail forming from the 
tidally disrupted star However, now the disk can form much 
closer to the black hole because the innermost circular orbit 
is at a smaller radius. Also, the higher rotation speed of the 
hole prevents rapid accretion of high angular momentum ma- 
terial. Thus the amount of nuclear matter remaining outside of 
the black hole is now much larger: about 30% of the neutron 
star mass. Within approximately 5 ms of disruption, a large 
fraction of that material (about half) forms an accretion disk 
(Fig. |5] center) with stable density and angular momentum 
profiles. 

The evolution of the resulting accretion disk is then simi- 
lar to what is expected for lower mass black holes. The disk 
expands slightly, and heats up to T ^ 5-6 McV, but does 



TABLE II: Parameters of the final configuration. i\/disk is defined 
as the mass available outside the black hole 5 ms after disruption. 
aBn/AfBH is the final dimensionless spin of the black hole, p™'"' is 
the maximum baryon density once the disk settles to an equilibrium 
configuration 20 ms after disruption. T is the average temperature 
and H is the semi-thickness of the disk at that time. 



Name 




A/bh 


max / g \ 

Po 


r(MeV) 


H_ 

r 




Q7S5 


< 0.004 


0.67 




Not Applicable 




86 


Q7S7 


0.06 


0.80 


2 X 10'" 


2 


0.3 


63 


Q7S9 


0.28 


0.92 


3 X 10'' 


6 


0.3 


39 


Q5S5 


0.06 


0.71 


8 X 10" 


2 


0.3 


106 



not evolve much otherwise. The disk is thick {H/R ^ 0.3), 
and an order of magnitude denser than for obh A^bh = 0.7. 
From the accretion rate at the end of the simulation, we de- 
duce an expected lifetime of about 75 ms, although this value 
would certainly be modified by the inclusion of magnetic ef- 
fects. Another important feature of this high-spin case is that 
a few percent of the total mass appears to be either unbound 
or on orbits such that it will not interact with the accretion 
disk within its expected lifetime. As the outermost part of the 
tidal tail is poorly resolved in our simulations (the grid points 
are focused in the region close to the black hole where the 
accretion disk forms), predictions about the future of material 
ejected far away from the black hole are not very reliable. The 
mere presence of such material, however, indicates that mas- 
sive ejecta might occur in BHNS mergers if the mass and spin 
of the black holes are high. By contrast, in previous simula- 
tions for lower BH masses or spins, all of the material in the 
tidal tail appeared unequivocally bound to the black hole. 

The large differences in the amount of matter available 
outside the black hole as we modify the black hole spin 
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FIG. 4: Matter distribution for simulation Q7S7. Left: Disruption of the neutron star. Material that is not immediately accreted onto the black 
hole (about 8% of theNS mass) forms a thin tidal tail. Density contours; po = (10^**, 10^^, 10^^, 10", 10^'')g/cm^. Center: Disk formation, 
about 8ms later. A low density accretion disk begins to form about 50 km away from the black hole, but rapid mass accretion from the tidal tail 
causes the disk profile to vary significantly over time. The maximum density is 3 x lO'^^ g/crci? . Density contours: po = (10^^, 10^^°) g/cm^. 
Right: About 20ms after disruption, a stable accretion disk has formed with low maximum density pmax ~ 2 x 10^" g/cm'^ . About 2.5% of the 
initial mass of the neutron star remains in the disk at that time, but accretion onto the black hole would destroy the disk within approximately 
20 ms. Density contour: po = 10*' g/cw? . 




FIG. 5: Matter distribution for simulation Q7S9. Left: Disruption of the neutron star, very similar to the Q7S7 case (Fig.|4l(, but with more 
material ejected into the tidal tail. Density contours: po ~ (10^^, 10^"^, 10^^, lO'^^, lO'^") g/cm"'. Center: Disk formation 6.5ms later. More 
than 25% of the fluid material remains outside the black hole. A massive accretion disk rapidly forms, containing more than 10% of the matter. 
Density contours: po ~ (10^^, 10^^, 10^") . Right: 20ms after disruption, a massive disk remains. It is 10 times denser than in the 

lower spin case, and has been in a nearly stable configuration for more than 10 ms. Density contours: po = (10^^ , lO^'') g/cm'\ 




can be seen in Fig. |6] The figure shows the difference be- 
tween a merger without disruption (Q7S5), with disruption 
followed by the slow formation of a low density, short-lived 
disk (Q7S7), and with disruption followed by the forma- 
tion of a heavy, mostly stable accretion disk (Q7S9). For 
A^bh/^j^ns = 7, we should thus expect strong qualitative dif- 
ferences in the behavior of BHNS mergers as we vary binary 
parameters such as the black hole spin or the NS equation of 
state, with a significant fraction of the parameter space now 
leading to direct merger of the binary without the formation of 
an accretion disk. From Fig.|6] we can also easily evaluate the 
influence of the black hole mass. The lower mass ratio simu- 
lation {q = 5), which has a black hole spin abh A^bh = 0.5, 
behaves approximately as the g = 7, aBH/.^^BH = 0.7 case. 



although the radius of the accretion disk around the less mas- 
sive black hole is naturally smaller, and the maximum baryon 
density much higher In Section |Vl we discuss in more detail 
how the qualitative behavior of BHNS mergers is affected by 
both the mass and spin of the black hole, and how this relates 
to their ability to power short gamma-ray bursts. 

The results of the mergers are summarized in Table [III In 
addition to the disk mass and properties discussed above. Ta- 
ble HI] lists the final spin of the BH and the velocity kick im- 
parted to the hole. As expected, the relative increase in the 
final black hole spin decreases as the initial spin of the black 
hole increases. Starting from a spin aBH/-^^BH = 0.5 leads 
to a final BH spin of 0.67, while an original spin of 0.9 only 
rises to 0.92 after merger For low spins, the velocity kicks 
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FIG. 6: Baryon mass available outside the black hole as a function of FIG. 7: Gravitational wave signal for a mass ratio q — 7, when vary- 
time. The mass is normalized by the initial mass of the neutron star. ing the black hole spin between aBn/A/BH = 0.5 - 0.9. The time 
1 ^ _ I I ^ I and phase of the waveforms are matched at the peak of the signal. 




-2 2 4 6 

t(ms) 

are close to binary black hole predictions ll46ll47ll . but as we 
go to higher spins, the magnitude of the kick decreases. This 
could be expected, as most of the contribution to the kick in bi- 
nary black hole mergers comes from the momentum radiated 
as gravitational waves just before merger, exactly the part of 
the signal that is not emitted by BHNS systems when the star 
is tidally disrupted before reaching the ISCO. The same effect 
was observed by Kyutoku et al. [70 at lower mass ratios. 



B. Gravitational wave signal 

The gravitational wave signal emitted as two compact ob- 
jects spiral in and merge contains information on the orbital 
parameters of the binary, the masses and spins of the compact 
objects, and, in the presence of a neutron star, on the equation 
of state of matter above nuclear density. In BHNS binaries, 
as in BH-BH systems, these parameters leave a mark on the 
waveform as the binary spirals in. But on top of the informa- 
tion from the orbital evolution of the binary, BHNS systems 
also show strong variations in the qualitative features of the 
gravitational wave signal at the time of merger. This is visible 
on Fig. I2I which shows the gravitational strain measured for 
various binaries with q ^ 7 but different spins, and Fig. [8] 
which shows the same signals in the frequency domain. The 
low spin simulation Q7S5 has a spectrum qualitatively simi- 
lar to that of a BH-BH binary : the power slowly decreases 
with increasing frequency as the binary spirals in, then peaks 
at the time of merger (~ 1 kHz), and finally falls off expo- 
nentially as the remnant black hole rings down. At higher 
spins (aBH/^j^BH > 0.7), the neutron star disrupts and we no 
longer observe a peak in the gravitational wave spectrum. The 
high-frequency signal now depends on the details of the tidal 
disruption of the star. Lower spins lead to less disruption, with 



FIG. 8: Power spectrum of the gravitational wave signal for a mass 
ratio 5 = 7, while varying the black hole spin between obh/A/bh = 
0.5 - 0.9. The cutoff at low frequency is only an effect of the finite 
length of the evolution. hcB is the effective amplitude of the gravita- 
tional wave signal, as defined in Eq. (41) of 



— Q7S7 




f(Hz) 



most of the mass still falling rapidly into the black hole, while 
higher spins lead to more significant disruption, and a more 
homogeneous spread of the neutron star material around the 
black hole. Accordingly, high spin spectra are cut off at lower 
frequency, but fall off more smoothly. Similar qualitative dif- 
ferences are found by Kyutoku et al. 17] and Etienne et al. jst). 
In fact. Fig. [8] is strikingly similar to Fig. 18 of which 
shows the spectra of the gravitational wave signal for BHNS 
binaries with q = 3 and spins asH /A/bh = —0.5, 0, 0.5. Ob- 
servationally, the effects of tidal disruption on the waveform 
are slightly easier to measure in the case of a high mass ratio, 
as the merger occurs at frequencies about a factor of 2 smaller 
than for g = 3 (/ ~ 1 - 1.5 kHz), but this still remains a 
challenging frequency range for Advanced LIGO. 

Figs.|9]and[T0]show the gravitational wave signal of the two 
o-bb/Mbh = 0.5 simulations in the time and frequency do- 
mains. As expected from the above discussion, the lower mass 
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FIG. 9: Gravitational wave signal for the low spin simulations 
(aBH/A/BH ~ 0.5) for mass ratios q = (5, 7). The time and phase 
of the waveforms are matched at the peak of the signal. 
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FIG. 10: Power spectrum of the gravitational wave signal for the low 
spin simulations (aBn/A/BH = 0.5) for mass ratios q = (5, 7). 
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ratio has a behavior similar to the q = 7, aBn/A/BH = 0.7 
case: both configurations have similar tidal disruption history, 
with a few percent of the mass being ejected in a tidal tail 
and the rest of the neutron star merging quickly with the black 
hole. The signal in the lower mass case (q = 5) is simply at a 
slightly higher frequency. 

It is worth noting that similar differences can also be due to 
changes in the radius of the neutron star ifyl- fioll . In order to 
extract information from that high-frequency signal, it is thus 
important to already have good estimates of the parameters of 
the binary to break that degeneracy. Given accurate gravita- 
tional wave templates, this can be obtained from the signal at 
lower frequency, during the orbital evolution. 



V. LIMITS ON DISK FORMATION AND GAMMA-RAY 
BURSTS 

All the simulations presented in this paper use fairly large 
stars, with i?NS — 14.4 km for A/ns = 1-4 Af©. They corre- 



spond to particularly favorable cases for disk formation. They 
thus allow us to deduce limits on the minimum spin required 
to form significant accretion disks, at least if we assume that 
real neutron stars have radii i?NS < 14.4 km (as predicted by 
Hebeler et al. ll45ll ). Numerical simulations of BHNS using 
nearly the same nuclear equation of state have already been 
carried at lower mass ratios by Etienne et al. as well as in 
our previous paper |@]. From these results, we notice that the 
behavior observed here in the cases q = 5, obh/A^bh = 0.5 
and q = 7, asH/A/BH = 0.7 is very similar to the results 
obtained for q ~ 3, obh/A^bh = 0.0 iHI^: the amount 
of matter remaining outside the black hole at late times is 
A/out/A^NS ^ 0.05 - 0.06. However, these three cases are 
not exactly equivalent, of course. The higher mass ratio cor- 
responds to a larger, lower density disk, which generally ap- 
pears less favorable to the generation of SGRBs. Nonethe- 
less, all the cases lead to a stable accretion disk of relatively 
low mass. Similarly, we have three cases for which nearly 
the entire star is quickly accreted onto the black hole, with a 
small but measurable amount of material ejected in a tidal tail 
A/out/A/NS ^ 0.005 - 0.008: the g = 7, obh/A^bh = 0.5 
case presented here, as well as the q ~ 5, obh/A/bh = 0.0 
and q ~ 3, obh/A/bh = —0.5 cases of Etienne et al. isj]. 
We thus have two curves in the space of black hole masses 
and spins giving us cases where the amount of matter avail- 
able at late time for a i?Ns = 14.4 km, A/ns = 1.4A'f0 star 
is - O.OIA/q and - O.OSA/q respectively (see Fig. [IB- Be- 
low the lowest curve of Fig. [TT] no significant disruption of 
the neutron star occurs before merger. Tidal effects are likely 
to affect the gravitational wave signal slightly, but the forma- 
tion of an accretion disk is not possible. Above the highest 
curve, massive disks are formed for this particular stellar ra- 
dius. However, all that can be said is that massive disks are 
possible: the more compact the neutron star, the more these 
curves would be displaced toward higher black hole spins. Fi- 
nally, in between the two curves, we form a low mass disk of 
fairly low density or just a tidal tail containing a few percent 
of the initial neutron star mass and slowly falling back onto 
the central black hole. 

In order for a BHNS binary to generate a short gamma-ray 
burst, a necessary condition would thus be to lie in the upper 
region of Fig.[TT] This is not, however, sufficient. More com- 
pact stars make it more difficult to form accretion disks(see 
e.g. ll48ll for a model based on earlier, lower mass ratio simu- 
lations), and the same is true if the spin of the black hole and 
the orbital angular momentum are not aligned, as shown in 
Foucart et al. |@] for q = 3 mass ratios. In i^, the effect of 
spin misalignment was fairly limited, but we were considering 
configurations for which even a nonspinning black hole led to 
the formation of a significant accretion disk. Nevertheless, the 
final disk mass appeared to depend on the aligned component 
of the black hole spin. There is no guarantee that this will still 
be true for higher mass ratio, or for nearly extremal spins; fur- 
ther numerical studies in a general relativistic framework will 
be necessary to determine the effect of spin misalignment in 
BHNS binaries for g ^ 7 - 10. But if the same pattern holds 
for q — 7 and q = 3, misalignment angles ^ 40° would be 
enough to have a BHNS with spin obh/A/bh = 0.9 behave 
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FIG. 1 1 : Simulations with similar mass remaining outside the black 
hole for various BH masses and spins. All simulations correspond 
to A'/ns ~ IAMq, _Rns = 14.4km (although a rescaling of A/ns, 
i?NS and A/bh would lead to the same fraction of the stellar mass re- 
maining at late time). They thus correspond to optimistic predictions 
for the formation of accretion disks. 
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as an aligned configuration with obhA/bh = 0.7. As mis- 
alignment angles of ^ 40° could be fairly common lITill . this 
would affect the probability that BHNS mergers form massive 
accretion disks. 

We should also note that BHNS binaries which do not form 
an accretion disk could still emit a detectable electromagnetic 
signal before the disruption of the neutron star, for example in 
the form of resonant shattering precursor flares ll49ll . 



VI. CONCLUSION 

Population synthesis models indicate that existing simu- 
lations of BHNS mergers in general relativity have proba- 
bly studied black holes with masses below the most common 
astrophysical values. Binaries with lower mass black holes 
could occur, and their study has already offered valuable in- 
sights into the influence of various binary parameters such as 
BH spin, NS equation of state and mass ratio on the emit- 
ted gravitational wave signal and the disk formation process. 
However, studies of the behavior of BHNS binaries for more 
massive black holes are necessary if we want to obtain more 
accurate predictions for what is likely to be a significant por- 
tion of the available BHNS parameter space. 

In this paper, we present the first simulations of BHNS bi- 
naries with A/bh = IOA/q. We show that, as opposed to what 
was found for lower mass black holes, the tidal disruption of 
the neutron star and the formation of an accretion disk are no 
longer the norm. We only observe significant accretion disks 
for high black hole spins, ubh/Mbh > 0.7, and this result is 



obtained using a fairly large neutron star (i?NS ~ 14.4 km), 
thus providing a likely lower bound on the spin required to 
form accretion disks in such systems. 

As spin-orbit misalignment also inhibits disk formation |@], 
it appears that the formation of massive accretion disks will 
be a common result of BHNS mergers only if large black hole 
spins aligned with the orbital angular momentum are frequent, 
or if black holes in BHNS binaries are less massive than ex- 
pected. This is of course of particular importance when deter- 
mining whether BHNS mergers can be the progenitors of short 
gamma-ray bursts. BHNS systems with A/bh ^ IOA/0 and 
o-bb/Mbu < 0.7 appear to be unlikely candidates to power 
SGRBs. In Sec.|V] we used existing simulations of BHNS bi- 
naries for mass ratios g = 3 - 7 to obtain approximate lower 
bounds for the magnitude of the black hole spin required to 
obtain massive disks. 

BHNS systems with Mbh > 7A/ns show another dif- 
ference with respect to previously studied binaries; for high 
black hole spins obhA^bh ^ 0.9, a few percent of the NS 
is either unbound or weakly bound. This could lead to the 
ejection of neutron-rich material in the neighborhood of the 
binary. However, study of this effect will require the use of 
numerical methods that resolve the evolution of the tidal tail 
better than our current code. 

We showed that we can extract gravitational waves with a 
cumulative phase eiTor Acf) < 0.2 rad over about five orbits, 
without applying any time or phase shift to the signal. When 
matching the time and phase of the signal at its peak, our ac- 
curacy is similar to that reported by Kyutoku et al. |01 for 
lower mass black holes. The gravitational wave signal shows 
that as the spin of the black hole decreases, the spectrum goes 
through three regimes. At high spins, the spectrum is cut off 
at low frequency as the star is disrupted and the remaining 
material forms a nearly axisymmetric accretion disk. For in- 
termediate spins, the spectrum is mostly flat, extending to fre- 
quencies ^ 2 kHz at which disruption leads to the formation 
of a thin tidal tail. At low spins, the spectrum is very similar to 
the signal from BH-BH binaries, where tidal disruption does 
not occur. 

The behavior of BHNS binaries for higher mass ratios 
{q > 7) or smaller stellar radii should be even more sim- 
ilar to BH-BH binaries, with small corrections for the tidal 
distortion of the star. Whether any such binary can form an 
accretion disk will depend on how mergers occur for quasi- 
extremal spins, and how frequent high-spin black holes are in 
BHNS binaries. For high-spin systems, the relative alignment 
of the black hole spin and the orbital angular momentum is 
also likely to play an important role. In an earlier study of 
BHNS binaries with misaligned spins we found that the 
qualitative effect of misalignment was relatively modest, but 
so were the effects of the spin itself, as even a nonspinning BH 
led to the formation of an accretion disk. At high mass ratios, 
it is likely that tidal tail and disk formation will only be pos- 
sible for small spin-orbit misalignments, as is also predicted 
by Newtonian simulations |19J. Additional simulations will, 
however, be necessary to obtain more accurate predictions for 
such systems. 
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